
********************************************************************************
*																			   *
*						Hidden Discount Regression Models					   *
*																			   *
********************************************************************************

/* 
 Note: To replicate must install zoib command "ssc install zoib". 
 Also must install spost13 commands from 
 https://jslsoc.sitehost.iu.edu/web_spost13/sp13_install.htm
*/

// Load data (ensure cd is set to "Data & Syntax" relative pathway)

cd "/Users/kpeter58/UNC Charlotte Dropbox/Kevin Petersen/Post-Doc/Active Discount Study/Data & Syntax"
use "Guideline_Data_CLEAN_FINAL.dta"

********************************************************************************
*						  Format and Label Variables						   *
********************************************************************************

// Create mean centered and fully standardized versions of numeric variables 
// Also create year sequence variable in case we want to look at time trends

sum Age_r
gen Age_0 = Age_r - r(mean)
egen Age_z = std(Age_r)

sum Low
gen Low_0 = Low - r(mean)
egen Low_z = std(Low)

sum Total_r
gen Total_0 = Total_r - r(mean)
egen Total_z = std(Total_r)

sum Total_r2
gen Total_02 = Total_r2 - r(mean)
egen Total_z2 = std(Total_r2)

sum N_Conv_r
gen N_Conv_0 = N_Conv_r - r(mean)
egen N_Conv_z = std(N_Conv_r)

********************************************************************************
*						  1. Main Effects Zoib Model 					       *
********************************************************************************

// Run model

zoib SSP_r2 i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult 

est store mod1

// Save out results for tables in R

putexcel set Results_FINAL.xlsx, sheet(Results) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

// MEMs for each panel, saving out for tables in R

margins, dydx(*) expression(predict(pr0)) atmeans nose

putexcel set Results_FINAL.xlsx, sheet(Pr0_MEM) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

margins, dydx(*) expression(predict(prcond)) atmeans nose

putexcel set Results_FINAL.xlsx, sheet(Prcond_MEM) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

margins, dydx(*) expression(predict(pr1)) atmeans nose

putexcel set Results_FINAL.xlsx, sheet(Pr1_MEM) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

// Predicted SSP for race, saving out for tables in R

mtable, at(Race_rc = (1 2 3)) atmeans statistics(all) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph1) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

********************************************************************************
*					      2. Race by Gender ZOIB Model						   *
********************************************************************************

zoib SSP_r2 i.Race_rc##i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc##i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc##i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc##i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult 

est store mod2

// Predicted SSP for race by gender, saving out for figures in R

mtable, at(Race_rc = (1 2 3) Sex_r = (1 2)) atmeans statistics(all) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph1) modify
	matrix A = r(table)'
	putexcel A9 = matrix(A), names
	
// Second differences table

qui mlincom (3 - 4) - (1 - 2), stat(est se z p) add rowname("Sex by race diff:Black - White")
	matrix A = (r(est), r(se), r(z), r(p))
qui mlincom (5 - 6) - (1 - 2), stat(est se z p) add rowname("Sex by race diff:Hispanic - White")
	matrix B = (r(est), r(se), r(z), r(p))
qui mlincom (3 - 4) - (5 - 6), stat(est se z p) add rowname("Sex by race diff:Black - Hispanic")
	matrix C = (r(est), r(se), r(z), r(p))
	
mlincom, title("Tests of Second Differences for Race and Gender") ///
	stat(est se z p) twidth(15)

// Combining matrices and saving out

matrix D = A\B\C
matrix colnames D = est se z p
matrix rownames D = Black-White Hispanic-White Black-Hispanic

putexcel set Second_Diffs_FINAL.xlsx, sheet(Race & Sex) modify
	putexcel A1 = matrix(D), names
	
********************************************************************************
* 				       3. Race by Gender and Age ZOIB Model                    *
********************************************************************************

zoib SSP_r2 c.Age_r##c.Age_r##i.Race_rc##i.Sex_r i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(c.Age_r##c.Age_r##i.Race_rc##i.Sex_r i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(c.Age_r##c.Age_r##i.Race_rc##i.Sex_r i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(c.Age_r##c.Age_r##i.Race_rc##i.Sex_r i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult

est store mod3

// Predicted SSP for race by gender and age, saving out for figures in R

mtable, at(Race_rc = (1 2 3) Sex_r = (1 2) Age_r = (18 21 30 40 50)) atmeans statistics(all) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph2) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names
	
// Third differences table

qui mlincom ((3 - 15) - (4 - 16)) - ((1 - 13) - (2 - 14)), stat(est se z p) add rowname("Age (18-30) by race/sex:Bm vs Bf - Wm vs Wf")
	matrix A = (r(est), r(se), r(z), r(p))
qui mlincom ((5 - 17) - (6 - 18)) - ((1 - 13) - (2 - 14)), stat(est se z p) add rowname("Age (18-30) by race/sex:Hm vs Hf - Wm vs Wf")
	matrix B = (r(est), r(se), r(z), r(p))
qui mlincom ((3 - 15) - (4 - 16)) - ((5 - 17) - (6 - 18)), stat(est se z p) add rowname("Age (18-30) by race/sex:Bm vs Bf - Hm vs Hf")
	matrix C = (r(est), r(se), r(z), r(p))

qui mlincom ((15 - 27) - (16 - 28)) - ((13 - 25) - (14 - 26)), stat(est se z p) add rowname("Age (30-50) by race/sex:Bm vs Bf - Wm vs Wf")
	matrix D = (r(est), r(se), r(z), r(p))
qui mlincom ((17 - 29) - (18 - 30)) - ((13 - 25) - (14 - 26)), stat(est se z p) add rowname("Age (30-50) by race/sex:Hm vs Hf - Wm vs Wf")
	matrix E = (r(est), r(se), r(z), r(p))
qui mlincom ((15 - 27) - (16 - 28)) - ((17 - 29) - (18 - 30)), stat(est se z p) add rowname("Age (30-50) by race/sex:Bm vs Bf - Hm vs Hf")
	matrix F = (r(est), r(se), r(z), r(p))
	
mlincom, title("Tests of Third Differences for Race by Gender and Age") ///
	stat(est se z p) twidth(15)
	
// Combining matrices and saving out

matrix G = A\B\C\D\E\F
matrix colnames G = est se z p
matrix rownames G = BmBfWmWf18 HmHfWmWf18 BmBfHmHf18 BmBfWmWf30 HmHfWmWf30 BmBfHmHf30

putexcel set Second_Diffs_FINAL.xlsx, sheet(Race Sex & Age) modify
	putexcel A1 = matrix(G), names
	
********************************************************************************
*					 4. Race by Mode of Conviction ZOIB Model 				   *
********************************************************************************

zoib SSP_r2 i.Race_rc##i.Disposition_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc##i.Disposition_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc##i.Disposition_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc##i.Disposition_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult 

est store mod4

// Predicted SSP for race by mode of conviction, saving out for figures in R

mtable, at(Race_rc = (1 2 3) Disposition_rc = (1 2 3 4)) atmeans statistics(all) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph3) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names
	
// Second differences table

qui mlincom (6 - 2) - (5 - 1), stat(est se z p) add rowname("Race nbp vs bp:B - W")
	matrix A = (r(est), r(se), r(z), r(p))
qui mlincom (10 - 2) - (9 - 1), stat(est se z p) add rowname("Race nbp vs bp:H - W")
	matrix B = (r(est), r(se), r(z), r(p))
qui mlincom (6 - 10) - (5 - 9), stat(est se z p) add rowname("Race nbp vs bp:B - H")
	matrix C = (r(est), r(se), r(z), r(p))
	
qui mlincom (7 - 3) - (5 - 1), stat(est se z p) add rowname("Race pna vs bp:B - W")
	matrix D = (r(est), r(se), r(z), r(p))
qui mlincom (11 - 3) - (9 - 1), stat(est se z p) add rowname("Race pna vs bp:H - W")
	matrix E = (r(est), r(se), r(z), r(p))
qui mlincom (7 - 11) - (5 - 9), stat(est se z p) add rowname("Race pna vs bp:B - H")
	matrix F = (r(est), r(se), r(z), r(p))
	
qui mlincom (8 - 4) - (5 - 1), stat(est se z p) add rowname("Race t vs bp:B - W")
	matrix G = (r(est), r(se), r(z), r(p))
qui mlincom (12 - 4) - (9 - 1), stat(est se z p) add rowname("Race t vs bp:H - W")
	matrix H = (r(est), r(se), r(z), r(p))
qui mlincom (8 - 12) - (5 - 9), stat(est se z p) add rowname("Race t vs bp:B - H")
	matrix I = (r(est), r(se), r(z), r(p))
	
mlincom, title("Tests of Second Differences for Race and Mode of Conviction") ///
	stat(est se z p) twidth(15)
	
// Combining matrices and saving out

matrix J = A\B\C\D\E\F\G\H\I
matrix colnames J = est se z p
matrix rownames J = nbpbpBW nbpbpHW nbpbpBH pnabpBW pnabpHW pnabpBH tbpBW tbpHW tbpBH

putexcel set Second_Diffs_FINAL.xlsx, sheet(Race & Dispo) modify
	putexcel A1 = matrix(J), names
		
********************************************************************************
*		     5. Race by Gender and Mode of Conviction ZOIB Model	     	   *
********************************************************************************

zoib SSP_r2 i.Race_rc##i.Sex_r##i.Disposition_rc Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc##i.Sex_r##i.Disposition_rc Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc##i.Sex_r##i.Disposition_rc Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc##i.Sex_r##i.Disposition_rc Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult

est store mod5 

// Predicted SSP for race by gender and mode of conviction, saving out for figures in R

mtable, at(Race_rc = (1 2 3) Sex_r = (1 2) Disposition_rc = (1 2 3 4)) atmeans statistics(all) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph3) modify
	matrix A = r(table)'
	putexcel A10 = matrix(A), names
	
// Third difference table

qui mlincom ((10 - 14) - (2 - 6)) - ((9 - 13) - (1 - 5)), stat(est se z p) add rowname("Sex/race/nbp vs bp:B - W")
	matrix A = (r(est), r(se), r(z), r(p))
qui mlincom ((18 - 22) - (2 - 6)) - ((17 - 21) - (1 - 5)), stat(est se z p) add rowname("Sex/race/nbp vs bp:H - W")
	matrix B = (r(est), r(se), r(z), r(p))
qui mlincom ((10 - 14)-(18 - 22)) - ((9 - 13)-(17 - 21)), stat(est se z p) add rowname("Sex/race/nbp vs bp:B - H")
	matrix C = (r(est), r(se), r(z), r(p))
	
qui mlincom ((11 - 15) - (3 - 7)) - ((9 - 13) - (1 - 5)), stat(est se z p) add rowname("Sex/race/pna vs bp:B - W")
	matrix D = (r(est), r(se), r(z), r(p))
qui mlincom ((19 - 23) - (3 - 7)) - ((17 - 21) - (1 - 5)), stat(est se z p) add rowname("Sex/race/pna vs bp:H - W")
	matrix E = (r(est), r(se), r(z), r(p))
qui mlincom ((11 - 15)-(19 - 23)) - ((9 - 13)-(17 - 21)), stat(est se z p) add rowname("Sex/race/pna vs bp:B - H")
	matrix F = (r(est), r(se), r(z), r(p))
	
qui mlincom ((12 - 16) - (4 - 8)) - ((9 - 13) - (1 - 5)), stat(est se z p) add rowname("Sex/race/t vs bp:B - W")
	matrix G = (r(est), r(se), r(z), r(p))
qui mlincom ((20 - 24) - (4 - 8)) - ((17 - 21) - (1 - 5)), stat(est se z p) add rowname("Sex/race/t vs bp:H - W")
	matrix H = (r(est), r(se), r(z), r(p))
qui mlincom ((12 - 16)-(20 - 24)) - ((9 - 13)-(17 - 21)), stat(est se z p) add rowname("Sex/race/t vs bp:B - H")
	matrix I = (r(est), r(se), r(z), r(p))
	
mlincom, title("Tests of Second Differences for Race by Sex and Mode of Conviction") ///
	stat(est se z p) twidth(15)
	
// Combining matrices and saving out

matrix J = A\B\C\D\E\F\G\H\I
matrix colnames J = est se z p
matrix rownames J = nbpbpBW nbpbpHW nbpbpBH pnabpBW pnabpHW pnabpBH tbpBW tbpHW tbpBH

putexcel set Second_Diffs_FINAL.xlsx, sheet(Race Sex & Dispo) modify
	putexcel A1 = matrix(J), names

********************************************************************************
*			6. Linear Models for Total and Active Sentence Length			   *
********************************************************************************

gen Total_ln = log(Total_r2)

regress Total_ln i.Race_rc i.Sex_r Age_0 c.Age_0#c.Age_0 i.Disposition_rc /* 
				*/ i.Prior_r i.Relationship_r i.Violations_r i.MaxSeverity_r Low_0 /*
				*/ i.Person_r i.Property_r i.Drug_r N_Conv_0 i.Weapon_r i.Minimum_r /*
				*/ i.Representation_rc i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r /*
				*/ i.Sentence_Year_r i.Jurisdiction, vce(robust)
				
// Saving out for tables in R
				
putexcel set Sensitivity_Mods_FINAL.xlsx, sheet(Log TS Length) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

// Predicted total sentence and active sentence using estimates from Figure 2

margins Race_rc, atmeans expression(exp(predict(xb))*exp((`e(rmse)'^2)/2)) post

display 76.31634 - (76.31634*0.764)
display 78.57414 - (78.57414*0.736)
display 72.22364 - (72.22364 *0.673)

putexcel set Predictions_FINAL.xlsx, sheet(Graph5) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

// Active sentence model with 95% CIs
// Adjusting the total sentence for each defendant by the predicted SSP
	
gen Active_pred = Total_r2 - (Total_r2*0.764) if Race_rc == 1
	replace Active_pred = Total_r2 - (Total_r2*0.736) if Race_rc == 2
	replace Active_pred = Total_r2 - (Total_r2*0.673) if Race_rc == 3
	
gen Active_predln = log(Active_pred + 1)
	
regress Active_predln i.Race_rc i.Sex_r Age_0 c.Age_0#c.Age_0 i.Disposition_rc /* 
				*/ i.Prior_r i.Relationship_r i.Violations_r i.MaxSeverity_r Low_0 /*
				*/ i.Person_r i.Property_r i.Drug_r N_Conv_0 i.Weapon_r i.Minimum_r /*
				*/ i.Representation_rc i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r /*
				*/ i.Sentence_Year_r i.Jurisdiction, vce(robust)
				
// The predictions from this model almost perfectly match crude predictions above
// Saving them out to create figure with error bars in R
				
margins Race_rc, atmeans expression(exp(predict(xb))*exp((`e(rmse)'^2)/2)) post

putexcel set Predictions_FINAL.xlsx, sheet(Graph5) modify
	matrix A = r(table)'
	putexcel A6 = matrix(A), names
	
********************************************************************************
*						   7. Race by Year ZOIB Model		       			   *
********************************************************************************

zoib SSP_r2 i.Race_rc##i.Sentence_Year_r i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.Disposition_rc i.Person_r i.Property_r i.Drug_r /*
			  */ i.MaxSeverity_r Low_z N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc##i.Sentence_Year_r i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.Disposition_rc i.Person_r i.Property_r i.Drug_r /*
		   */ i.MaxSeverity_r Low_z N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc##i.Sentence_Year_r i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.Disposition_rc i.Person_r i.Property_r i.Drug_r /*
			*/ i.MaxSeverity_r Low_z N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Jurisdiction) /*
*/ phivar(i.Race_rc##i.Sentence_Year_r i.Sex_r Age_z c.Age_z#c.Age_z i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.Disposition_rc i.Person_r i.Property_r i.Drug_r /*
	   */ i.MaxSeverity_r Low_z N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Jurisdiction) /* 
*/ robust difficult 

est store mod6

// MEMs for race by year, saving out for figures in R

forval i = 1/24 {
	mtable, dydx(Race_rc) at(Sentence_Year_r = `i') atmeans statistics(all)
	
	putexcel set Years_w_CI_FINAL.xlsx, sheet(Year`i') modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names
}

// Second differences table (2 = year 2000)

mtable, at(Race_rc = (1 2 3) Sentence_Year_r = (2 24)) atmeans statistics(all) post

qui mlincom (4 - 2) - (3 - 1), stat(est se z p) add rowname("Race 2023 vs 2020:B - W")
	matrix A = (r(est), r(se), r(z), r(p))
qui mlincom (6 - 2) - (5 - 1), stat(est se z p) add rowname("Race 2023 vs 2020:H - W")
	matrix B = (r(est), r(se), r(z), r(p))
qui mlincom (4 - 6) - (3 - 5), stat(est se z p) add rowname("Race 2023 vs 2020:B - H")
	matrix C = (r(est), r(se), r(z), r(p))
	
mlincom, title("Tests of Second Differences for Race by Year") ///
	stat(est se z p) twidth(15)
	
// Combining matrices and saving out

matrix D = A\B\C
matrix colnames D = est se z p
matrix rownames D = BW20232000 HW20232000 BH20232000

putexcel set Second_Diffs_FINAL.xlsx, sheet(Race & Year) modify
	putexcel A1 = matrix(D), names
	
********************************************************************************
*				  8. Sensitivity Analysis with Adult Sample		       		   *
********************************************************************************

frame copy default Adult
frame change Adult
drop if Age_r < 23

zoib SSP_r2 i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult


putexcel set Sensitivity_Mods_FINAL.xlsx, sheet(Adult) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names
	
********************************************************************************
*			 9. Sensitivity Model Including Credit and Home Detention		   *
********************************************************************************

frame change default

replace SSP_r = 0 if SSP_r < 0

zoib SSP_r i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.Race_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult 


putexcel set Sensitivity_Mods_FINAL.xlsx, sheet(Credit & Detention) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names

********************************************************************************
*			   10. Sensitivity Model w/ Race and Ethnicity Combos	    	   *
********************************************************************************

recode RaceEth_rc (6 = 2) (2 = 4) (4 = 6)
label define RaceEth 1 "White_NH" 2 "White_H" 3 "Black_NH" 4 "Black_H" 5 "Other_NH" 6 "Other_H"
label values RaceEth_rc RaceEth

zoib SSP_r2 i.RaceEth_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			  */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			  */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			  */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction, /*
*/ oneinflate(i.RaceEth_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
		   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
		   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
		   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ zeroinflate(i.RaceEth_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
			*/ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
			*/ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
			*/ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /*
*/ phivar(i.RaceEth_rc i.Sex_r Age_z c.Age_z#c.Age_z i.Disposition_rc i.Prior_r i.Relationship_r /* 
	   */ i.Violations_r i.MaxSeverity_r Low_z i.Person_r i.Property_r i.Drug_r /*
	   */ N_Conv_z i.Weapon_r i.Minimum_r i.Representation_rc /*
	   */ i.Victim_r i.Injury_r i.Impact_r i.Vulnerability_r i.Sentence_Year_r i.Jurisdiction) /* 
*/ robust difficult

putexcel set Sensitivity_Mods_FINAL.xlsx, sheet(Race & Ethnicity) modify
	matrix A = r(table)'
	putexcel A1 = matrix(A), names




